Bird beak's elaboration and exploration

All the custom code and custom functions called in this document can be downloaded and installed in R as the standalone working package beer (Beak Elaboration and Exploration in R). Although beer is intentionally designed to be portable and shareable, we advice workers to use it with moderation and tailor it's consumption to their specific research needs. Some analyses in beer do take some time to run so we advice to not drive or operate heavy machinery while using it. Don't use it while driving. The following code snippets are used to illustrate the implementation of the functions. See for running the actual analyses.

Elaboration and exploration analyses

## Loading required package: dispRity

## 
## Attaching package: 'beer'

## The following object is masked from 'package:datasets':
## 
##     trees

Extracting the variance covariance matrices from the MCMCglmm mini-chains

To analyse the exploration and elaboration aspects of both clades and tips, we first selected a random 1000 covariance matrices from the combined MCMCglmm posterior variance-covariance matrices.

set.seed(42)
## Selecting 1000 random covariance matrices from the MCMCglmm
covar_matrices <- get.covar(combined_results, n = 1000)[1:4]

We will use these covariance matrices as the basis of the phylogenetic main axes of variation. These axes correspond to the main axes of variation in the \(N_DIMENSIONS\) trait space at different levels:

• The whole bird phylogeny level ("phylogeny level");
• And the named node levels ("clade levels").

These major axes where calculated as the longest distance within the 95% confidence interval hyper-ellipse from the variance-covariance matrix (i.e. the major axes of the 95% CI ellipse of matrix). In practice, the coordinates of these axes end points were calculated directly from the cross product between the 95% unit hypersphere scaled by the variance-covariance matrix's squared eigen values and the variance-covariance matrix's eigen vectors.

Major axes calculation

For \(O_{n}\), the unit hypersphere matrix of n dimensions and of radius being composed of the two identity matrices \(I_{n}\) and \(-I_{n}\) so that:

\[\begin{equation} O_{n} = \begin{pmatrix} 1 & 0 & \cdots & 0 \\ 0 & 1 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & 1 \\ -1 & 0 & \cdots & 0 \\ 0 & -1 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & -1 \\ \end{pmatrix} \end{equation}\]

(in other words \(O_{n}\) is the matrix representing each edge point of the unit hyper sphere of \(n\) dimensions and of radius \(1\));

And \(O'_{n}\) the scaled matrix hyper sphere to the right confidence interval size using the \(\chi^2\) distribution:

\[\begin{equation} O'_{n} = O_{n} \sqrt{\chi^2(0.95)} \end{equation}\]

For the variance-covariance matrix \(VCV_{n}\) of \(n\) dimensions

\[\begin{equation} VCV_{n} = \begin{pmatrix} \sigma(a) & \sigma(a,b) & \cdots & \sigma(a,n) \\ \sigma(a,b) & \sigma(b) & \cdots & \sigma(b,n) \\ \vdots & \vdots & \ddots & \vdots \\ \sigma(n,a) & \sigma(n,b) & \cdots & \sigma(n) \\ \end{pmatrix} \end{equation}\]

and the eigenvectors v and the eigenvalues \(\lambda\) satisfying the following eigen decomposition:

\[\begin{equation} VCV_{n} \bold{v} = \lambda \bold{v} \end{equation}\]

We can get \(M_{n}\) matrix containing all the edges' coordinates of the 0.95 CI hypersphere from \(VCV_{n}\) using the transposition of the cross product between the eigenvectors v and the product of the scale 0.95 CI unit sphere \(O'_{n}\) and the eigenvalues \(\lambda\):

\[\begin{equation} M_{n} = [(O'_{n}\sqrt{\lambda}) \times \bold{v}]^{\text{T}} \end{equation}\]

Where \(M_{1,n}\) is the major axes of the 0.95 CI hyper ellipse fitting the variance-covariance matrix, \(M_{2,n}\) the second second axes (minor axes in 2D), etc...

Finally, the centred the matrix \(M_{1,m}\) on the centroid of each clade.

The detailed procedure was adapted from 李哲源's post on Stack Overflow and implemented in beer::get.axes (specific procedure on lines 40-73).

## Get the centre of each clade
group_centres <- list(
    "gulls"      = colMeans(demo_data[demo_data$clade == "gulls", c(1:4)]),
    "plovers"    = colMeans(demo_data[demo_data$clade == "plovers", c(1:4)]),
    "sandpipers" = colMeans(demo_data[demo_data$clade == "sandpipers", c(1:4)]),
    "phylogeny"  = colMeans(demo_data[, c(1:4)]))

## Calculating the major axes
major_axes <- get.axes(covar_matrices, centre = group_centres)

One can visualise these ellipses and major axes in 2D using the beer::plot.axes and beer::plot.ellipses functions (based on Murdoch and Chow (2020)):

set.seed(51)
## Sub-selecting 5 random VCV matrices for readibility
sub_sample_matrices <- get.covar(combined_results, 5)[1:4]
sub_sample_major_axes <- get.axes(sub_sample_matrices, centre = "intercept")

## Plotting the ellipses
plot.ellipses(sub_sample_matrices, col = c("grey",colour_vector),
              add = FALSE, transparent.scale = 0.7)
plot.axes(sub_sample_major_axes, col = c("grey",colour_vector),
              add = TRUE, transparent.scale = 0.7)
legend("topleft", legend = c("phylogeny", levels(demo_data$clade)), col = c("grey", colour_vector), lty = 1)

We can also visualise all the major axes in relation to the elements in the space.

## Plotting it
plot.space(demo_data,
           col = colour_vector,
           levels = demo_data$clade,
           xlab = "PC1 (90.5%)",
           ylab = "PC1 (6.86%)",
           xlim = c(-2, 2),
           ylim = c(-1, 1))
points(do.call(rbind,group_centres)[,c(1,2)],
       pch = 13, col = c(colour_vector, "grey"))
plot.axes(major_axes["animal"], col = "grey",
          add = TRUE, transparent.scale = 0.01,
          centre = 0)
plot.axes(major_axes["animal:clade1"], col = colour_vector[1],
          add = TRUE, transparent.scale = 0.01,
          centre = colMeans(demo_data[demo_data$clade == "gulls", c(1,2)]))
plot.axes(major_axes["animal:clade2"], col = colour_vector[2],
          add = TRUE, transparent.scale = 0.01,
          centre = colMeans(demo_data[demo_data$clade == "plovers", c(1,2)]))
plot.axes(major_axes["animal:clade3"], col = colour_vector[3],
          add = TRUE, transparent.scale = 0.01,
          centre = colMeans(demo_data[demo_data$clade == "sandpipers", c(1,2)]))
legend("topleft", legend = c("phylogeny", levels(demo_data$clade)), col = c("grey", colour_vector), lty = c(1,1,1,1), pch = c(NA, 19, 19, 19))

We can then compare these axes either between each other or relative to the tips positions in the trait space using elaboration and exploration scores:

Elaboration and exploration scores calculations

From these major axes, we can use linear algebra in \(n\) dimensions to estimate clades or tips' exploration and elaboration scores. These scores are a direct interpretation of Endler et al. (2005)'s elaboration and innovation (here exploration) patterns (Fig 1. in Endler et al. (2005)): Given a major axis of evolutionary "trajectory" (in the sense of evolutionary trend, not an evolutionary direction), on can calculate the elaboration score as the projection value of any vector defined by a clade or a tip on the major axis and the exploration (or "innovation" in Endler et al. (2005)) as that vector's rejection.

Figure: in a space with 5 elements: A, B, C, D, E (in grey); where D and E represent the major axes (e.g. the major phylogenetic effect axes from the MCMCglmm), we can rotate and rescale each element so that D and E become the unit vector of length 1 (the black letters D and E) and get the exploration elaboration scores for either the elements (e.g. element A's projection in blue and its rejection in orange) or from any other axes (e.g. the major phylogenetic effect axes from the MCMCglmm from a specific clade).

Measuring projection and rejection

For any vectors \(\vec{a}\) and \(\vec{b}\), defined either as one set or a pair of sets of coordinates in \(n\) dimensions:

\[\begin{equation} \vec{a} = \begin{bmatrix} x \\ y \\ \cdots \\ n \\ \end{bmatrix} \end{equation}\] \[\begin{equation} \vec{a} = \begin{bmatrix} x_{1} & x_{2} \\ y_{1} & y_{2} \\ \cdots & \cdots \\ n_{1} & n_{2} \\ \end{bmatrix} \end{equation}\]

We can calculate \(\vec{a_{1}}\), the orthogonal projection of \(\vec{a}\) onto \(\vec{b}\) using:

\[\begin{equation} \vec{a_{1}} = \frac{\vec{a} \cdot \vec{b}}{\|\vec{b}\|} \end{equation}\]

With \(\|\vec{b}\| = \sqrt{\vec{b} \cdot \vec{b}}\) being the norm of \(\vec{b}\). And \(\vec{a_{2}}\), the rejection of \(\vec{a}\) onto \(\vec{b}\):

\[\begin{equation} \vec{a_{2}} = \vec{a} - \vec{a_{1}} \end{equation}\]

Generalisation of projection onto any vector in a set space

Using this, we can calculate the projection and rejection for any element within a trait space \(TS_{m,n}\):

\[\begin{equation} TS_{m,n} = \begin{bmatrix} x_{1} & x_{2} & \cdots & x_{m} \\ y_{1} & y_{2} & \cdots & y_{m} \\ \vdots & \vdots & \ddots & \vdots \\ n_{1} & n_{2} & \cdots & n_{m} \\ \end{bmatrix} \end{equation}\]

And any base vector \(\vec{B}\) defined as:

\[\begin{equation} B = \begin{bmatrix} x_{1} & x_{2}\\ y_{1} & y_{2}\\ \vdots & \vdots \\ n_{1} & n_{2} \\ \end{bmatrix} \end{equation}\]

By using the linear transformation \(f_{\vec{B}}\) of the trait space \(TS\) moving \(\vec{B}\) onto \(TS\)'s first axes unit vector \(\vec{\hat{\imath}}\):

\[\begin{equation} f_{\vec{B}}(TS) = \left( \frac{TS - [Bx_{1}, By_{1}, \cdots, Bn_{1}]^{\text{T}}}{\|\vec{B}\|} \right) \cdot R_{\vec{B}} \end{equation}\]

With R_{} being the rotation matrix of the vector \(\vec{B}\) onto \(\vec{\hat{\imath}}\):

\[\begin{equation} R_{\vec{B}} = I_{\vec{B}} - \vec{B}\vec{B}^\text{T} - \vec{\hat{\imath}}\vec{\hat{\imath}}^\text{T} + [\vec{B} \vec{\hat{\imath}}] \begin{bmatrix} cos(\theta) & -sin(\theta)\\ sin(\theta) & cos(\theta)\\ \end{bmatrix} [\vec{B} \vec{\hat{\imath}}]^\text{T} \end{equation}\]

Where \(\theta\) is:

\[\begin{equation} \theta = acos \left(\frac{\vec{B} \cdot \vec{\hat{\imath}}}{\|\vec{B}\| \cdot \|\vec{\hat{\imath}}\|} \right) \end{equation}\]

Or \(\theta = acos (B_x)\) since both \(\|\vec{B}\|\) and \(\|\vec{\hat{\imath}}\|\) are equal to 1 and \(\|\vec{\hat{\imath}}\|\) is the unit vector on the first axis.

Algorithm for calculating the projection/rejection of any element in a defined space

In practice we followed this procedure and applied a modification of this implementation (see Aguilera and Pérez-Aguila (2004) for the formal generalisation of this algorithm) using the following algorithm:

1. In the trait space, define \(\vec{B}\) as the base vector (typically \(\vec{B}\) is defined as the pair of coordinates set from the major axis described above).
2. Centre the trait space on the origin of \(\vec{B}\) so that the first set of coordinates of \(\vec{B}\) are 0.
3. Scale the trait space to the norm of \(\vec{B}\) so that the norm of \(\vec{B}\) is now 1.
4. Rotate the trait space using the rotation matrix \(R_{\vec{B}}\) to satisfy the linear transformation \(\vec{B} \arrow \vec{\hat{\imath}}\) (with \(\vec{\hat{\imath}}\) being the first unit vector of the trait space - typically the x axes unit vector).
5. Project/reject every element in the trait on \(\vec{B}\) (that is now \(\vec{\hat{\imath}}\)). In practice, the first coordinate (x) of each element is now it's projection onto \(\vec{B}\).

Following the ways to: 1. estimate the variance-covariance matrices of the phylogenetic effects in trait space (02-MCMCglmm_mini_chains.Rmd); 2. extracting the major axes of variation of these VCV matrices for the different levels in the phylogeny; 3. and calculating the elaboration and exploration scores for elements in the trait space based on these major axes;

We can measure these scores for each clade and each species separately.

Elaboration and exploration by clade

To run the analyses per clade we compared the major axes for each clade to the major axes for the overall phylogenetic effect (e.g. plovers' VCV vs. all birds' VCV, etc.).

To do so, for randomly selected posterior variance-covariance matrix we defined \(\vec{a}\) as the major axes for the clade of interest and projected it onto \(\vec{B}\) as the major overall phylogenetic effect (see above). For each pair of \(\vec{a}\) and \(\vec{B}\), we translated the origin of \(\vec{a}\) onto the origin of \(\vec{B}\) and measured the projection (elaboration), rejection (exploration) and angle (degree of exploration).

## Get all the angles and stuff
group_results <- analyses.group(major_axes, base = "animal")
names(group_results) <- levels(demo_data$clade)

The results are reported as boxplots for each group indicating which group explores/elaborates more (and to which degree - angle) relatively to each other.

## Plot the blob wise metrics
plot.analyses.group(group_results, col = colour_vector)

Elaboration and exploration by species

Similarly to the analyses described above, we ran the analyses per species but here we directly compared each species position in the space to the same vector \(\vec{B}\) described above and reported the results per groups.

## Global analyses
tip_results_global <- analyses.tip(data = demo_data[, c(1:4)], axes = major_axes$animal)

And here are the values plotted by group:

## Plot results
plot.analyses.tip(tip_results_global, col = colour_vector, group = split(rownames(demo_data), f = demo_data$clade))

Note that here only the tips positions projected on the global phylogenetic axes.

We also looked at the exploration elaboration profiles within each clades by replacing for each clade, the projection of the elements onto \(\vec{B}\) to the projection onto their respective clade major axes \(\vec{a}\).

TODO!

References

Aguilera, Antonio, and Ricardo Pérez-Aguila. 2004. “General N-Dimensional Rotations.” Václav Skala-UNION Agency.

Endler, John A, David A Westcott, Joah R Madden, and Tim Robson. 2005. “Animal Visual Systems and the Evolution of Color Patterns: Sensory Processing Illuminates Signal Evolution.” Evolution 59 (8). Wiley Online Library: 1795–1818.

Murdoch, Duncan, and E. D. Chow. 2020. Ellipse: Functions for Drawing Ellipses and Ellipse-Like Confidence Regions. https://CRAN.R-project.org/package=ellipse.